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Abstract 

In many applications, there is a desire to determine if the dynamics of interest 
are chaotic or not. Since positive Lyapunov exponents are a signature for chaos, they 
are often used to determine this. Reliable estimates of Lyapunov exponents should 
demonstrate evidence of convergence; but literature abounds in which this evidence 
lacks. This paper presents two maps through which it highlights the importance 
of providing evidence of convergence of Lyapunov exponent estimates. The results 
suggest cautious conclusions when confronted with real data. Moreover, the maps 
are interesting in their own right. 
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1 Introduction 

Computations of Lyapunov exponents continue to play a significant role in the study 
of nonlinear systems. The computations are often numerical because most nonlinear 
systems prove analytically intractable. In the case of real systems, there may be no 
mathematical model that adequately describes the system; hence one is often confronted 
with a limited amount of data from which to estimate Lyapunov exponents. Finite ob- 
servations may afford computations of finite time Lyapunov exponents, but not global 
Lyapunov exponents. Unlike global Lyapunov exponents, finite time Lyapunov expo- 
nents depend on the initial conditions This limits their value in characterising the 
underlying system. 

Fortunately, [2] provides a theorem that guarantees convergence of finite time Lya- 
punov exponents to the global Lyapunov exponents. The hypothesis of the theorem 
(also found in pQ) requires the underlying system to have an invariant measure. It 
behoves us, therefore, to satisfy ourselves that the system of interest has an invariant 
measure and that our estimates converge to the global Lyapunov exponent. Yet exam- 
ples abound in which Lyapunov exponent estimates are provided with neither of these 
two requirements confirmed (e.g. [3j[U[5l[6j[7l[8l[9]). Some have emphasised the need 
to provide confidence limits along with the estimates [lOl E] , but confidence limits are 
valuable only when there is evidence of convergence. 

This paper presents two maps through which it highlights the importance of demon- 
strating evidence of convergence of Lyapunov exponent estimates. One map has a pa- 
rameter that takes non-negative integer values. For all possible values of this parameter, 
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we can compute the invariant measure in closed form. The other map presents computa- 
tional challenges, despite its illusive simplicity. The two maps are presented in § [3j with 
computations of Lyapunov exponents. We conclude with a discussion and summary of 
the results in § [H Lyapunov Exponents are briefly discussed in the following section. 

2 Lyapunov Exponents 

Lyapunov exponents measure the average rate of separation of two trajectories that are 
initially infinitely close to each other. Consider a map 



where xt € R m - For an initial state xq, the dynamics of its small perturbation, Sxq, are 
governed by the linear propagator, M.(xq, N), which is a product of Jacobians so that 
M(xq, N) = DF(xn-x) ■ ■ ■ DF{xi)DF{xq). The finite time average separation/growth 
rate of two initially nearby trajectories is then given by 



Oseledec [2] provides a theorem that guarantees that the limit limjv->oo Ajv is unique. 
Denote this limit by A. If 5xq is a member of the right singular vectors of A4, then Atv 
is a finite time Lyapunov exponent and A is a global Lyapunov exponent. 

It has been noted that m > 1 corresponds to non-commuting matrices [HIH]. On 
the other hand, when m = 1, we take logs of positive scalars so that the ergodic theorem 
may be used to yield 



where p(x) is the invariant distribution of the dynamics. If a dynamical system settles 
onto an invariant distribution, then the effect of transients on dynamical invariants 
averages out. Moreover, equation ([1]) provides an alternative to computing Lyapunov 
exponents when m = 1; hence, simple bootstrap resampling techniques may be applied 
in a straight forward way to make Lyapunov exponent estimates, with no need for block 
resampling approaches suggested by Ziehmann et al. [11] . Given a set of observations, 
bootstrapping is accomplished by repeatedly sampling randomly from the data with 
replacement to estimate the statistic of interest [11]. This would give a distribution of 
the statistic of interest, which in this case is A^v- 

Since, for a given N, finite time Lyapunov exponents are functions of the initial 
states, it is useful to report values corresponding to a distribution of initial states to 
determine convergence. For an assessment of convergence of finite time Lyapunov ex- 
ponents, one can sample the initial states from the invariant distribution, p(x). In nu- 
merical computations, the invariant distribution is often not accessible in closed form. 
Nonetheless one can iterate forward an initial distribution po(x) so that after N itera- 
tions pn{x) is an estimate of the invariant distribution. One would then sample initial 
states from pn{x) to estimate Ajy's. The aim here is to sample the initial states accord- 
ing to the invariant measure, if it exists. As a consequence of Oseledec's theorem [2], the 
distribution of the Atv's will converge to a delta distribution centred at A as N — > oo. 

When one is faced with real data, the approach of the previous paragraph cannot be 
used unless there is a reliable mathematical model of the system. Nonetheless, one may 
use bootstrap approaches suggested in [11] for various values of N to assess convergence 
in the distribution of the Aat's. 



x t+ i = F(x t ) 
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Figure 1: (Left) Graph of the Infinitely piece- wise linear map given in equation ([2]). 



3 The Two Maps 



In this section, two novel maps with contrasting behaviours are considered. Analytic 
computations are complemented with numerical results and the importance of estab- 
lishing convergence is highlighted. 

3.1 Infinitely Piece-wise Linear Map 

Let us consider the infinitely piece-wise linear map 



defined for fixed k G {0, 1, . . . , oo}. A graph of the map corresponding to k = 4 is shown 
in figure [TJ The Lyapunov exponents for this map may be sought either analytically or 
numerically. For a given value of A;, an exact knowledge of the invariant measure of this 
map is sufficient for one to determine the corresponding global Lyapunov exponent. 

3.1.1 Invariant Densities 

Denote the invariant density corresponding to a given value of k by p( k \x) and the 
associated invariant measure H by For k = 0, the invariant distribution of this map 
is the uniform distribution U[0, 1]. That is 




l-i 



(2) 



P 




l, ze[0,i], 

0, otherwise. 
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To see this, we note that the invariant density must satisfy the Perron-Frobenius equa- 
tion [12j [13] , whence 



,(o) 



i=i 



It is then evident that the uniform distribution satisfies the above equation identically. 
When k > 0, the following propositions hold and their proofs are given in appendix [Al 

Proposition 1 For k = 1, the invariant density of the piecewise linear map is 

1 -e[0,*), 



3' 



3' 



Proposition 2 For any k > 2, the invariant density of the infinitely piecewise linear 
map is 

:l - X G [0, p-) , 



/ 2 fc +2 



p (fc) (x) 



2 

3' 

2 
3' 



X ^ [ 2 fc ' 2) ' 
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3.1.2 Global Lyapunov Exponents 

The Lyapunov exponents of the map are given by 

-1 



p (fc) (x)log 2 0' fc (x)dx, 



where the derivative is taken only at points of continuity, and is thus a weak derivative. 
We then state the following proposition whose proof is given in appendix [Bl 

Proposition 3 For any k G {0, 1,2,...}, the global Lyapunov exponent of the infinitely 
piecewise linear map is 

= 2. 



3.1.3 Numerical Results 

Numerical estimates of global Lyapunov exponents may be obtained in either of two 
ways. One approach is to compute finite time Lyapunov exponents along trajectories 
from different initial conditions. The distribution of finite time Lyapunov exponents 
should converge to a delta distribution as the length of the trajectories increases. For 
the infinitely piecewise linear map in equation ([2]) with k = 4, we considered trajectories 
of length N = 2 n starting from 2 10 initial conditions. We varied n from to 20 and 
the resulting distributions are shown in figure [2j Notice that with every doubling of 
the length of trajectories, the distribution of Lyapunov exponents converges to a delta 
distribution that is centred at 2. 

Alternatively, we can plot finite time Lyapunov exponents corresponding to each ini- 
tial condition as a function of trajectory length, N. These are also plotted in figure [2] on 
the right hand side. All the lines corresponding to 2 4 initial conditions clearly converge 
to a value around 2. 
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Figure 2: (Left) Distributions of finite time Lyapunov exponents for the infinitely piecewise map 



in ([2]). The different distributions correspond to different time lengths of trajectories 
according to the colour bar. (Right) Finite time Lyapunov exponents from 2 4 as a 
function of time, N = 2 n . The different lines correspond to different initial conditions. 



Finally, we could compute the Lyapunov exponents by first estimating the invariant 
density. The invariant density may be estimated by evolving forward the standard uni- 
form distribution ?7[0, 1] under the map 4>k( x )- Again we considered k = 4 and evolved 
forward M points sampled from U[0, 1]. A time series of distributions corresponding to 
M = 2 16 is shown in figure [3j The distribution appears to have converged to the invari- 
ant distribution after N = 2 n , when n = 7. Lyapunov exponent estimates corresponding 
to M = 2 13 (left) and M = 2 16 (right) are shown in figure HI where the solid line cor- 
responds to the mean and the dash-dotted lines are bootstrap uncertainty estimates. 
Notice that the value of A^ fc ^ = 2 is generally within the confidence intervals for n > 4. 
The confidence intervals on the left are relatively big. This is due to M = 2 13 being too 
small to provide a good estimate of the invariant distribution. As may be noted from 
the right hand graphs on figure [21 convergence was obtained with more than N = 2 15 
points. It is, therefore, not surprising that when M = 2 16 the confidence intervals are 
much smaller as is the case on the right hand graph of the same figure. 

3.2 Paradoxical Map 

The map considered here is 




2(1 -x), xe (1/2,1). 



xe (0,1/2), 



(3) 



Its graph is given in figure [5] on the left panel. It has an unstable fixed point at the 
origin. The stability of the origin is determined using 
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Figure 3: (Left) Evolution of 2 16 points sampled from U[0, 1] under the infinitely piecewise map 
in ([2]) with the number of iterations given by N — 2" . The histograms were obtained 
by using 2 4 equally spaced bins on the unit interval. 



Clearly, tp (x) > 1 for any x > 0. However, when x is close to zero then <p (x) is close 
to 1. Hence the origin is a weak repeller as pointed out in [12J. Indeed it is this property 
that makes the dynamics of this map paradoxical as was found in [12\ 113] on a similar 
map that differs only by the linear part. If ^{x) € (0, 1/2) for all j = 0, 1, . . . , n — 1, 
then 



1 



nx 



where <p°(x) = x. On the hand, if ^{x) £ (1/2, 1) for all j = 0, 1, . . . , n — 1, then 



tp n (x) 



1 - (-2) n 



(-2) 



n-l. 



It is, therefore, evident that with probability one trajectories with initial conditions sam- 
pled according to the Lebesgue measure on [0, 1] do not converge to the origin. However, 
it can be shown that iterates of distributions that are initially uniform converge to a 
delta distribution centred at the origin. These two contrasting behaviours of trajectories 
and densities are paradoxical and are due to the weak repeller at the origin. 

Numerically, the invariant density of this map may be determined by evolving for- 
ward in time the uniform distribution U [0,1]. Iterating forward an initial ensemble of 
2 15 points sampled from U[0, 1] yields the graphs shown in figure [5] on the right panel. 
From the graphs, (p N (U[0, 1]) appears to have nearly converged to the invariant density 
when N = 2 24 . If that is indeed the case, it follows that one can compute estimates of 
the global Lyapunov exponents using points sampled from ip N (U[0, 1]) when N = 2 25 
via equation (pQ). Unfortunately, the dynamics are trapped into very long periodic orbits 
after about N = 2 25 iterations, the most dominant period being of the order of 2 29 . 

Distributions of finite time Lyapunov exponents for varying N are shown in figure [6] 
on the right panel. The initial conditions of trajectories used to estimate each Ajv were 
sampled from ip N (U[0, 1]). Whereas it seems the distribution of Lyapunov exponents has 
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Figure 4: Lyapunov exponent estimates from a distribution obtained by evolving forward M = 
2 13 (left) and M = 2 16 (right) points sampled from U[0, 1] under the infinitely piece- 
wise linear map when k = 4. Each distribution is thus (U[0, 1]), where N = 2™. In 
each graph, the solid line corresponds to the mean estimate and the dash-dotted lines 
correspond to bootstrap uncertainty estimates. 

converged when N = 2 25 , the distribution converged to is not a delta distribution. In 
fact, if iterations are continued further, the distribution ultimately becomes tri-modal, 
each mode corresponding to each periodic orbit that the dynamics settle onto. If one 
had only a finite sample of data without a knowledge of the data-generating process, 
as is typical in practical situations, they could mistakenly attribute the bell shape to 
uncertainty. 

Finally, graphs of Ajv ; computed along trajectories starting from different initial 
conditions, versus log N indicated that there is indeed no convergence (see figure [6] on 
the left panel). Each line on the graph corresponds to a number of initial conditions, M. 
This raises suspicion over Lyapunov exponent estimates obtained from a single initial 
condition. An example of where this was done is [4], albeit using real data. 

4 Discussion 

The aim of this paper was to highlight the importance of demonstrating convergence 
of finite time Lyapunov exponents estimates of the global Lyapunov exponent. Two 
contrasting maps were thus presented, one of which its dynamical properties are known 
accurately, and the other not. The former map has a parameter that takes non- negative 
integer values. As the sole parameter is varied across all non-negative integer values, 
the global Lyapunov exponent remains fixed at A = 2. At a sample of parameter 
values, numerical computations confirmed a fast convergence of distributions of finite 
time Lyapunov exponents to a delta distribution. The centre of the delta distributions 
coincided with the analytically found value of 2. 

The other map was more challenging. Detailed insights into its dynamics were sought 
numerically. Numerical distributions of finite time Lyapunov exponent estimates did not 
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Figure 5: (Left) Graph of the map <p(x) given in equation (J3)). (Right) Evolution of 2 15 ini- 
tially uniformly distributed ensemble in [0, 1] under the map (p(x). The bins for these 
histograms were logarithmically spaced to the base 2. 

converge to a delta distribution; even with the use of trajectories as long as 2 25 . If the 
length of trajectories is increased, the dynamics ultimately fall onto very long periodic 
orbits due to finite machine precision. The lack of convergence of distributions of finite 
time Lyapunov exponents to a delta distribution should alarm us against reporting the 
mode as an estimate of the global Lyapunov exponent. Indeed any arising confidence 
limits would equally be nonsensical. In fact, following analytic considerations of a similar 
map in |12t [T3] would indicate that the map has no invariant distribution. Hence, for a 
randomly selected initial condition, the dynamics are transient with probability one. 

The failure of distributions of finite time Lyapunov exponents to converge to a delta 
distribution is a caution against reporting estimates of Lyapunov exponents and con- 
fidence limits. In numerical computations, convergence of the distributions to a delta 
distribution is achievable if the dynamics settle onto an invariant distribution before 
being trapped onto a numerical periodic orbit. A higher machine precision may be nec- 
essary to ensure this. When using data from a real nonlinear system, stationarity has 
to be established first. Convergence of distributions of finite time Lyapunov estimates 
should then be verified by progressively increasing the length of trajectories from which 
the estimates are made. Block resampling approaches suggested in (TTJ may help provide 
the distributions that are used to assess convergence to a delta distribution. 
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Figure 6: (left)Graphs of finite time Lyapunov exponents versus trajectory length. 

(right) Distributions of Lyapunov Exponents. Each distribution is obtained 
by computing finite time Lyapunov exponents along trajectories whose initial 
conditions are sampled from U[0, 1]. 
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A Proofs of Propositions 1 and 2 

Proof of Proposition 1: When k = 1, we note that the invariant density is piecewise 
constant with 

r ci, i£ii = [0,|) , 

( c 2 , x £ A 2 = [|,l] , 

where c\ and c 2 are constants. In order to determine the constants c\ and c 2 , we note 
that by definition an invariant measure satisfies the property // fc ) (A) = ^{^{A)} [JJ, 
where (f>^ n {A) is the set of all the points that are mapped onto A after n iterations of 
the map 4>k(x). It thus follows that 

^){A 1 )= l -^)(A l ) + ^{A 2 ). (A.l) 

For a probability measure, it is also true that 

^\A x ) + ^ l \A 2 ) = l. (A.2) 

Solving ([AH) and (|A~2"1) yields 
Hence c\ = 4/3 and c 2 = 2/3. 

Proof of Proposition 2: Let us now consider the general case of k > 2. In this 
case, the invariant density is constant on three intervals. That is 





f r (k) 

c l > 


x G A[ k) = 




p {k \x) = < 


M 


x G A (k) = 






(fe) 


x G A (k) = 


[*.!]■ 
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Again, we use the underlying invariant measure to compute the c - s. If we let fj,j 
H {k \Aj), then 



(k) 
Mi 



(fe) 

M2 



(fe) 
M 3 



1 

2 k V 



(k) . (Jfc) 

m! + m 2 



, (ft) 

+ M 3 > 



1 1 

2 ~ 2^ 



(fc) . (k) 
Ml + M 2 



(k) . (k) 
Mi + M 2 



(A.3) 
(A.4) 
(A.5) 



Notice that the sum of the first two equations yields the last equation. Therefore, we 

(k) 

need another equation in order to determine the fij s, which is 



(fc) . (k) . (k) , 
Mi + M 2 + M 3 = !■ 

Using any two of equations (|A.3[) , (|A.4p and (|A.5P and equation (|A.6P yields 



(A.6) 



„<*> _ + 1 „W 
1 ~ 3 • 2 fc ~ 1 ' 2 



jjfe-i 



3 • 2 



k-l 



(k) 

M 3 



From these it follows that 



( fc )_2* + 2 {k) 2 



2 ~3> c 3 



(fe) = 2 
3' 



B Proof of Proposition 3 

The proof of this proposition considers each of the cases k = 0, k = 1 and k > 2 
separately. When A; = 0, we get 



A (o) 



1=1 2' 1 
oo 

£ 



P (0) log 2 <Ao(x)da; 
log 2 (p' (x)dx 

log 2 </>o(z)d:z: 



00 „ 1 

2'~ 



1=1 

2. 
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Similarly, for k = 1 we get 
A« = 



f p ( - 1 \x)\og 2 (f> 1 (x)dx 
Jo 

o / log 2 0i(x)dx + - / log 2 1 (x)dx 
•5 Jo -3 Jl/2 



. oo 

4 x - i 

i=2 



= 2. 



For general k > 2, 

A« = / 1 p( fc )(x)log 2 ^(x)d 

JO 



cf ) / 2 log 2 0fc(x)dx + 4 fc) / log 2 ^(x)da; + 4 fc) / log 2 fc (x)dx 
7o ^ 



(2 fc + 2) 



*\ ; 2 i 

i=fc+l i=2 



(2 fe + 2) (A; + 2) 2 1 



+ = 2. 



In the above derivation, we used the identity 



Ei _ (m + 1) 
oi om-1 
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